package statistics;

import java.math.BigDecimal;
import edu.rit.numeric.Statistics;

public class D_Contrasts extends AbstractStatisticTests implements IStatisticalTest {
	
	static BigDecimal[] lnFactorialList = new BigDecimal[5000];
	static {
		lnFactorialList[0] = new BigDecimal("0");
		for(int i=1; i<5000; i++) {
			BigDecimal number = new BigDecimal( Math.log( i ) );
			lnFactorialList[i] = lnFactorialList[i-1].add(number);
		}
	}
	
	
	/*
	 * Heavy tailed approximation of Normal for Case-only analysis 
	 */
	public double TestRareCombinations ( int order, double lambda, int observed ) 
	{	
		double poisson_p =0;
		
		for(int k=observed; k<=observed+10; k++) {
			BigDecimal numerator = new BigDecimal( - lambda + k*Math.log(lambda) );
			double log_poisson_p = ( numerator.subtract( lnFactorialList[k] ) ).doubleValue(); //e^-\lambda * lambda^k / k!
			poisson_p += Math.exp(log_poisson_p);
		}
		
		return poisson_p;
	}
	
	
	/*
	 * (non-Javadoc)
	 * @see statistics.AbstractStatisticTests#TestCasesOnly(int, int[][][], int[][][], int)
	 */
	public double TestCasesOnly ( int order, int[][][] observedCases ) 
	{			
		double z_cases = TestCasesOnly_Zscore(order, observedCases);
		double	probability = Statistics.chiSquarePvalue(1, Math.pow(z_cases,2)) ; // / 2;
		return probability;
	}
	
	/*
	 * (non-Javadoc)
	 * @see statistics.AbstractStatisticTests#TestCasesVsControls(int, int[][][], int[][][], int, int[][][], int[][][], int)
	 */
	public double TestCasesVsControls ( int order, int[][][] observedCases, int[][][] observedControls ) 
	{ 
		double z_diff = TestCasesVsControls_Zscore(order, observedCases, observedControls);
		double probability = Statistics.chiSquarePvalue(1, Math.pow(z_diff,2)) ; // / 2;
		return probability;
	};
	
	
	
	public double TestCasesOnly_Zscore ( int order, int[][][] observedCases ) 
	{	
		int totalCases = observedCases[1][1][1] + observedCases[1][1][0] + observedCases[1][0][1] + observedCases[1][0][0] ;
		
		double P = (double) ( observedCases[1][1][1] + observedCases[1][1][0] ) / totalCases;		
		double Q = (double) ( observedCases[1][1][1] + observedCases[1][0][1] ) / totalCases;
		
		double expectedCoCarriers = P*Q*totalCases;
		double D_cases = (double) (observedCases[1][1][1] -  expectedCoCarriers) / totalCases;
		double sigma_cases  = Math.sqrt( P*(1-P)*Q*(1-Q) / totalCases );
		
		double z_cases = D_cases / sigma_cases;
		return z_cases;
	}
	
	
	public double TestCasesVsControls_Zscore ( int order, int[][][] observedCases, int[][][] observedControls ) 
	{ 
		int TOTAL_CASES = observedCases[1][1][1] + observedCases[1][1][0] + observedCases[1][0][1] + observedCases[1][0][0] ;
		int TOTAL_CONTROLS = observedControls[1][1][1] + observedControls[1][1][0] + observedControls[1][0][1] + observedControls[1][0][0] ;
		double MIN_CASE_FREQ = (double) 1 / (TOTAL_CASES + 1);
		double MIN_CONTROL_FREQ = (double) 1 / (TOTAL_CONTROLS + 1);
		
		//	Calculate case deviation
		double P = (double) ( observedCases[1][1][1] + observedCases[1][1][0] ) / TOTAL_CASES;		
		double Q = (double) ( observedCases[1][1][1] + observedCases[1][0][1] ) / TOTAL_CASES;
		P = Math.max(P, MIN_CASE_FREQ);
		Q = Math.max(Q, MIN_CASE_FREQ);
		double expectedCaseCarriers = P*Q*TOTAL_CASES;
		
		double D_cases = (double) (observedCases[1][1][1] -  expectedCaseCarriers) / TOTAL_CASES;	
		double sigma_cases  = Math.sqrt( P*(1-P)*Q*(1-Q) / TOTAL_CASES );

		//	Calculate control deviation
		double p = (double) ( observedControls[1][1][1] + observedControls[1][1][0] ) / TOTAL_CONTROLS;
		double q =  (double) ( observedControls[1][1][1] + observedControls[1][0][1] ) / TOTAL_CONTROLS;
		p = Math.max(p, MIN_CONTROL_FREQ);
		q = Math.max(q, MIN_CONTROL_FREQ);
		double expectedControlCarriers = p*q*TOTAL_CONTROLS;
		
		double D_controls = (double) (observedControls[1][1][1] -  expectedControlCarriers) / TOTAL_CONTROLS;
		double sigma_controls = Math.sqrt( p*(1-p)*q*(1-q) / TOTAL_CONTROLS );		
		
		//	LD-contrast
		double D_diff = D_cases - D_controls; 
		double sigma_diff = Math.sqrt( Math.pow(sigma_cases, 2)  + Math.pow(sigma_controls, 2) );
		double z_diff = D_diff / sigma_diff;	
		
		return z_diff;
	}
	
}
